In this script I’ll be looking at how to analyse the outputs from the MCMCglmm models to get the something about the elaboration/exploration stories.

I’m going to look at it using two approaches: * the blob-wise one where I’ll look at the differences between whole blobs (ellipses) in terms of elaboration/exploration (Robinson & Beckerman style); * and the tip-wise one where I’ll look at the elaboration/exploration score of the tips on the tree (cf. their blobs) relative to the whole phylogeny and to their respective blobs (clades).

Data

I’m going to focus on the data from Gavin and especially the models 5 and 6 that are:

## Loading the correct models
load("../Data/Processed/model_list.rda")
model_phylo1_clade3 <- model_list[[4]]
model_phylo3_clade3 <- model_list[[6]]

And here’s what the trait space looks like:

## Loading the data
load("../Data/Processed/morphdat.rda")
## Plotting it
colour_vector <- c("orange", "blue", "darkgreen")
plot(morphdat[, c(1,2)], pch = 19,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
legend("topleft", legend = levels(morphdat$clade), col = colour_vector, pch = 19)

Note that the PC% are from the whole PC in Gavin’s example.

Blob-wise approach

Here we’re basically comparing multidimensional ellipses (here 3D ones but I think we should push it to more dimensions!).

Phylo + clade model

First we can visualise some of these ellipses (100, randomly drawn from the posterior):

source("../Functions/plot.ellipses.R")
source("../Functions/get.covar.R")

## Selecting 100 random covariance matrices from the MCMCglmm
covar_matrices <- get.covar(model_phylo1_clade3,
                            n = 100,
                            simplify = FALSE)

## Plotting the results
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "none",
              col = c("grey",colour_vector),
              add = TRUE)

Because this is really messy, we can centre these ellipses on the groups ellipses average centres:

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "level",
              col = c("grey",colour_vector),
              add = TRUE)

We can then measure different stuff from these ellipses (e.g. going full statistics from Robinson & Beckerman).

However, here I’m just going to focus on the elaboration/exploration of these ellipses. For the three clades, I’m going to measure: * the angle of the clade’s ellipse to the general one (which will be their degree of exploration); * their scaled projection (where the ellipses all start from the point 0,0): showing how much the clade elaborates on the major phylo axis; * and their scaled rejection (how much they explore).

For these three metrics I expect: * the gulls and plovers to elaborate quiet a lot (angle ~90) compared to the sandpipers (angle ~90) * the gulls to explore the more followed by the plovers and then the sandpipers * the sandpipers to ellaborate the more followed by the gulls and then the plovers

Note that the three metrics are measured independently of the ellipse position in space and are measured in 3D.

source("../Functions/get.axes.R")
## Get all the phylo main axes
level_centred_major_axes <- get.axes(covar_matrices, centre = "level")
uncentered_major_axes <- get.axes(covar_matrices, centre = "none")

And this is what it looks like for the uncentred ones:

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "none",
              col = c("grey",colour_vector),
              add = TRUE)
plot.all.axes(uncentered_major_axes,
         col = c("grey",colour_vector),
         add = TRUE)

And for the centred ones:

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "level",
              col = c("grey",colour_vector),
              add = TRUE)
plot.all.axes(level_centred_major_axes,
         col = c("grey",colour_vector),
         add = TRUE)

Or without all the points n stuff

par(mfrow = c(1,2))
plot.all.axes(uncentered_major_axes,
         col = c("grey",colour_vector),
         add = FALSE, main = "Major axes (un-centred)")
plot.all.axes(uncentered_major_axes,
         col = c("grey",colour_vector),
         add = FALSE, main = "Major axes (level-centred)")

Phylo clade + clade model

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo3_clade3, n = 100,
              centre = "none", add = TRUE, col = c("black", colour_vector, colour_vector))

This looks much more messier (and also probably not scaled correctly?)…

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo3_clade3, n = 100,
              centre = "level", add = TRUE, col = c("black", colour_vector, colour_vector))

Tip-wise approach